#Last updated on 7th Feb, 2014

#TABLE 1#######################################################################################################################

rm(list=ls())
ls()

#Installing packages
install.packages("QuantPsyc")
install.packages("pwr")

#loading data
library(foreign)
a<-read.spss("Data.sav",to.data.frame=TRUE)
dim(a)#161 participants 334 variables

a<-a[-94,]#romoved case #94 with intellectual disabilities since it wasn't included in the whole process of analysis
dim(a)#160 participants 334 variables

#in the paper, verbal fluency was the average raw score from 5 tasks
a$verbal.fluency<-a$verbal.fluency/5
a$verbal.fluency

#in the paper, nonverbal fluency was the average raw score from 2 tasks
a$design.fluency<-a$design.fluency/2
a$design.fluency


#dividing participants into 3 subgroups (T:typically developing, L:low language functioning, S:specific language impairment) 
T<-subset(a,JCPPGroup=="T                       " )
dim(T)#88 participants for T group
L<-subset(a,JCPPGroup== "LLF                     " )
dim(L)#31 participants for L group
S<-subset(a,JCPPGroup=="S                       ")
dim(S)#41 participants for S group


#Firstly, calculating means & sd of three groups
#In the paper,18 variables were used out of 334 variables
#The mean of S(specific language impairment group)
S<-S[,c("age.test","bas.mat.t","bas.v.std","rs.scale","fs.scale","wc.rec.scale","wc.exp.scale","celf.age.eq","wmtbc.lr.raw","ooo.total","verbal.fluency","design.fluency","vimiVtot.err","vimiMtot.err","dkef.st.v.raw","dkef.st.nv.raw","TMswitch.cost","ied.tot.err")]
Smean<-round(sapply(S,mean,na.rm=TRUE),1)
Smean

#The sd of S(specific language impairment group)
Ssd<-round(sapply(S,sd,na.rm=TRUE),1)
Ssd

#The mean of L(low language functioning group)
L<-L[,c("age.test","bas.mat.t","bas.v.std","rs.scale","fs.scale","wc.rec.scale","wc.exp.scale","celf.age.eq","wmtbc.lr.raw","ooo.total","verbal.fluency","design.fluency","vimiVtot.err","vimiMtot.err","dkef.st.v.raw","dkef.st.nv.raw","TMswitch.cost","ied.tot.err")]
Lmean<-round(sapply(L,mean,na.rm=TRUE),1)
Lmean

#The sd of L(low language functioning group)
Lsd<-round(sapply(L,sd,na.rm=TRUE),1)
Lsd

#The mean of T(typically developing group)
T<-T[,c("age.test","bas.mat.t","bas.v.std","rs.scale","fs.scale","wc.rec.scale","wc.exp.scale","celf.age.eq","wmtbc.lr.raw","ooo.total","verbal.fluency","design.fluency","vimiVtot.err","vimiMtot.err","dkef.st.v.raw","dkef.st.nv.raw","TMswitch.cost","ied.tot.err")]
Tmean<-round(sapply(T,mean,na.rm=TRUE),1)
Tmean

#The sd of T(typically developing group)
Tsd<-round(sapply(T,sd,na.rm=TRUE),1)
Ssd

#TABLE 1(means and sds)
S_L_T<-cbind(Smean,Ssd,Lmean,Lsd,Tmean,Tsd)
S_L_T


colnames(S_L_T)<-c("SLI_mean","SLI_sd","LLF_mean","LLF_sd","Typical_mean","Typical_sd")
rownames(S_L_T)<-c("Age_months","BAS Matrices T-score","BAS verbal IQ","Recalling sentences","Formulated sentences","Word classes receptive","Word classes expressive","Language age_months","ELWM verbal","ELWM nonverbal","Fluency verbal","Fluency nonverbal","Inhibition verbal","Inhibition nonverbal","Planning verbal","Planning nonverbal","Switching verbal","Switching nonverbal")
S_L_T

#TABLE 1 (ranges)
S<-S[,c("age.test","bas.mat.t","bas.v.std","rs.scale","fs.scale","wc.rec.scale","wc.exp.scale","celf.age.eq","wmtbc.lr.raw","ooo.total","verbal.fluency","design.fluency","vimiVtot.err","vimiMtot.err","dkef.st.v.raw","dkef.st.nv.raw","TMswitch.cost","ied.tot.err")]
Srange<-round(sapply(S,range,na.rm=TRUE),1)
Srange

L<-L[,c("age.test","bas.mat.t","bas.v.std","rs.scale","fs.scale","wc.rec.scale","wc.exp.scale","celf.age.eq","wmtbc.lr.raw","ooo.total","verbal.fluency","design.fluency","vimiVtot.err","vimiMtot.err","dkef.st.v.raw","dkef.st.nv.raw","TMswitch.cost","ied.tot.err")]
Lrange<-round(sapply(L,range,na.rm=TRUE),1)
Lrange

T<-T[,c("age.test","bas.mat.t","bas.v.std","rs.scale","fs.scale","wc.rec.scale","wc.exp.scale","celf.age.eq","wmtbc.lr.raw","ooo.total","verbal.fluency","design.fluency","vimiVtot.err","vimiMtot.err","dkef.st.v.raw","dkef.st.nv.raw","TMswitch.cost","ied.tot.err")]
Trange<-round(sapply(T,range,na.rm=TRUE),1)
Trange





#TABLE 2##################################################################################################################################################################
# Hierarchical regression
# reg(odd number, e.g.,reg1,reg3,reg5...)is Step 1 regression using 2 predictive variables (age, nonverbal IQ).
# reg(even number, e.g., reg2, reg4, reg6...)is Step 2 regression using 4 predictive variables (age, nonverbal IQ, SLI-vs.-LLF group, SLE-vs.-typcial group).

attach(a)
#Step1 regression for ELWM verbal
reg1<-lm(wmtbc.lr.raw~age.test+bas.mat.t)
summary(reg1)

#Step2 regression for ELWM verbal
reg2<-lm(wmtbc.lr.raw~age.test+bas.mat.t+dummy2.sli.llf+dummy1.sli.typ)
summary(reg2)


library(QuantPsyc)
#converting coefficients to standardised beta values for Step 2 regression for ELWM verbal
#rounding standardised bata values to the second decimal place
b<-round(lm.beta(reg2),2)
p <- summary(reg2)$coefficients[-1,4]  

#adding asterisks to the p-values of coefficients
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
b<-paste(b, mystars, sep="")

#Step1 regression for ELWM nonverbal
reg3<-lm(ooo.total~age.test+bas.mat.t)
summary(reg3)

#Step2 regression for ELWM nonverbal
reg4<-lm(ooo.total~age.test+bas.mat.t+dummy2.sli.llf+dummy1.sli.typ)
summary(reg4)
c<-round(lm.beta(reg4),2)
p <- summary(reg4)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
c<-paste(c, mystars, sep="")

#Step1 regression for Fluency verbal
reg5<-lm(verbal.fluency~age.test+bas.mat.t)
summary(reg5)

#Step2 regression for Fluency verbal
reg6<-lm(verbal.fluency~age.test+bas.mat.t+dummy2.sli.llf+dummy1.sli.typ)
summary(reg6)
d<-round(lm.beta(reg6),2)
p <- summary(reg6)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
d<-paste(d, mystars, sep="")

#Step1 regression for Fluency nonverbal
reg7<-lm(design.fluency~age.test+bas.mat.t)
summary(reg7)

#Step2 regression for Fluency nonverbal
reg8<-lm(design.fluency~age.test+bas.mat.t+dummy2.sli.llf+dummy1.sli.typ)
summary(reg8)
e<-round(lm.beta(reg8),2)
p <- summary(reg8)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
e<-paste(e, mystars, sep="")

#Step1 regression for Inhibition verbal
reg9<-lm(vimiVtot.err~age.test+bas.mat.t)
summary(reg9)

#Step2 regression for Inhibition verbal
reg10<-lm(vimiVtot.err~age.test+bas.mat.t+dummy2.sli.llf+dummy1.sli.typ)
summary(reg10)
f<-round(lm.beta(reg10),2)
p <- summary(reg10)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
f<-paste(f, mystars, sep="")

#Step1 regression for Inhibition nonverbal
reg11<-lm(vimiMtot.err~age.test+bas.mat.t)
summary(reg11)

#Step2 regression for Inhibition nonverbal
reg12<-lm(vimiMtot.err~age.test+bas.mat.t+dummy2.sli.llf+dummy1.sli.typ)
summary(reg12)
g<-round(lm.beta(reg12),2)
p <- summary(reg12)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
g<-paste(g, mystars, sep="")

#Step1 regression for Planning verbal
reg13<-lm(dkef.st.v.raw~age.test+bas.mat.t)
summary(reg13)

#Step2 regression for Planning verbal
reg14<-lm(dkef.st.v.raw~age.test+bas.mat.t+dummy2.sli.llf+dummy1.sli.typ)
summary(reg14)
h<-round(lm.beta(reg14),2)
p <- summary(reg14)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
h<-paste(h, mystars, sep="")

#Step1 regression for Planning nonverbal
reg15<-lm(dkef.st.nv.raw~age.test+bas.mat.t)
summary(reg15)

#Step2 regression for Planning nonverbal
reg16<-lm(dkef.st.nv.raw~age.test+bas.mat.t+dummy2.sli.llf+dummy1.sli.typ)
summary(reg16)
i<-round(lm.beta(reg16),2)
p <- summary(reg16)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
i<-paste(i, mystars, sep="")

#Step1 regression for Switching verbal
reg17<-lm(log(TMswitch.cost)~age.test+bas.mat.t)
summary(reg17)

reg18<-lm(log(TMswitch.cost)~age.test+bas.mat.t+dummy2.sli.llf+dummy1.sli.typ)
summary(reg18)
j<-round(lm.beta(reg18),2)
p <- summary(reg18)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   "))))  
j<-paste(j, mystars, sep="")

#Step1 regression for Switching nonverbal
reg19<-lm(ied.tot.err~age.test+bas.mat.t)
summary(reg19)

#Step2 regression for Switching nonverbal
reg20<-lm(ied.tot.err~age.test+bas.mat.t+dummy2.sli.llf+dummy1.sli.typ)
summary(reg20)
k<-round(lm.beta(reg20),2)
p <- summary(reg20)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
k<-paste(k, mystars, sep="")

#a matrix of standardised beta values of Step 2 regression for 10 variables
Dmat<-rbind(b,c,d,e,f,g,h,i,j,k)
Dmat

#adding column 1 for "Total R^2 accounted for by the model" & column 6 for "��R^2 Step2"
Dmat<-data.frame(a=c(0,0,0,0,0,0,0,0,0,0),Dmat,a=c(0,0,0,0,0,0,0,0,0,0))
Dmat

#column 1: R^2 values from summaries of reg(even numbers), which are Step 2 regressions for each of 10 variables
Dmat[,1]<-c(round(summary(reg2)$r.squared,2),
round(summary(reg4)$r.squared,2),
round(summary(reg6)$r.squared,2),
round(summary(reg8)$r.squared,2),
round(summary(reg10)$r.squared,2),
round(summary(reg12)$r.squared,2),
round(summary(reg14)$r.squared,2),
round(summary(reg16)$r.squared,2),
round(summary(reg18)$r.squared,2),
round(summary(reg20)$r.squared,2))
Dmat

#column 2: the difference of R squared between Step 1 and Step 2 for each of 10 variables
anova(reg1,reg2)

#obtaining the p-value of anova test of reg1 and reg 2
p <- anova(reg1,reg2)[-1,6] 

#addingg asterisks to the p-values of anova test
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 

#adding asterisks to the difference of R squared between reg 1 and reg 2
p1<-paste(round((summary(reg2)$r.squared)-(summary(reg1)$r.squared),2), mystars, sep="")

#the followings show the same procedures (as shown above for the for reg 1 and reg 2 for the first variable) for obtaining the difference of R squared between Step 1 and Step 2 for each variable

p <- anova(reg3,reg4)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p2<-paste(round((summary(reg4)$r.squared)-(summary(reg3)$r.squared),2), mystars, sep="")


p <- anova(reg5,reg6)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p3<-paste(round((summary(reg6)$r.squared)-(summary(reg5)$r.squared),2), mystars, sep="")


p <- anova(reg7,reg8)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p4<-paste(round((summary(reg8)$r.squared)-(summary(reg9)$r.squared),2), mystars, sep="")


p <- anova(reg9,reg10)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p5<-paste(round((summary(reg10)$r.squared)-(summary(reg9)$r.squared),2), mystars, sep="")


p <- anova(reg11,reg12)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p6<-paste(round((summary(reg12)$r.squared)-(summary(reg11)$r.squared),2), mystars, sep="")


p <- anova(reg13,reg14)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p7<-paste(round((summary(reg14)$r.squared)-(summary(reg13)$r.squared),2), mystars, sep="")

p <- anova(reg15,reg16)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p8<-paste(round((summary(reg16)$r.squared)-(summary(reg15)$r.squared),2), mystars, sep="")


p <- anova(reg17,reg18)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p9<-paste(round((summary(reg18)$r.squared)-(summary(reg17)$r.squared),2), mystars, sep="")


p <- anova(reg19,reg20)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p10<-paste(round((summary(reg20)$r.squared)-(summary(reg19)$r.squared),2), mystars, sep="")


#column 6 shows that Step 2 explains more of the overall variance than Step 1 model
#��R^2 = R^2 of Step2 - R^2 of Step1
#making a matrix of the differences of R squared between Step 1 and Step 2 for 10 variables
Dmat[,6]<-c(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10)

#putting column/row names 
colnames(Dmat)<-c("Total R^2 accounted for by the model","age","NVIQ","SLI-vs.-LLF","SLI-vs.-typical","R^2 step2")
rownames(Dmat)<-c("ELWM verbal","ELWM nonverbal","Fluency verbal","Fluency nonverbal","Inhibition verbal","Inhibition nonverbal","Planning verbal","Planning nonverbal","Switching verbal","Switching nonverbal")

Dmat


#TABLE3############################################################################################################################################################################################################
# analyses in the same way as for TABLE2 but with one additional predictive variable, verbal IQ, in Step 1 and Step 2 regression analyses
# reg(odd number, e.g.,reg1,reg3,reg5...)is Step 1 regression using 3 predictive variables (age, nonverbal IQ, verbal IQ)
# reg(even number, e.g., reg2, reg4, reg6...)is Step 2 regression using 4 predictive variables (age, nonverbal IQ, verbal IQ, SLI-vs.-LLF group, SLE-vs.-typcial group)

# Hierarchical regression
attach(a)

#Step1 regression for ELWM verbal
reg1<-lm(wmtbc.lr.raw~age.test+bas.mat.t+bas.v.std)
summary(reg1)

#Step2 regression for ELWM verbal
reg2<-lm(wmtbc.lr.raw~age.test+bas.mat.t+bas.v.std+dummy2.sli.llf+dummy1.sli.typ)
summary(reg2)

install.packages("QuantPsyc")
library(QuantPsyc)

#standardised beta values of Step 2 regression for ELWM verbal
b<-round(lm.beta(reg2),2)

#adding asterisks to the p-vaules of coefficients of Step 2 regression results
p <- summary(reg2)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
b<-paste(b, mystars, sep="")

#Step1 regression for ELWM nonverbal
reg3<-lm(ooo.total~age.test+bas.mat.t+bas.v.std)
summary(reg3)
lm.beta(reg1)

#Step2 regression for ELWM nonverbal
reg4<-lm(ooo.total~age.test+bas.mat.t+bas.v.std+dummy2.sli.llf+dummy1.sli.typ)
summary(reg4)
c<-round(lm.beta(reg4),2)

#adding asterisks to the p-vaules of coefficients of Step 2 regression results
p <- summary(reg4)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
c<-paste(c, mystars, sep="")

#Step1 regression for Fluency verbal
reg5<-lm(verbal.fluency~age.test+bas.mat.t+bas.v.std)
summary(reg5)

#Step2 regression for Fluency verbal
reg6<-lm(verbal.fluency~age.test+bas.mat.t+bas.v.std+dummy2.sli.llf+dummy1.sli.typ)
summary(reg6)
d<-round(lm.beta(reg6),2)

#adding asterisks to the p-vaules of coefficients of Step 2 regression results
p <- summary(reg6)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
d<-paste(d, mystars, sep="")

#Step1 regression for Fluency nonverbal
reg7<-lm(design.fluency~age.test+bas.mat.t+bas.v.std)
summary(reg7)

#Step2 regression for Fluency nonverbal
reg8<-lm(design.fluency~age.test+bas.mat.t+bas.v.std+dummy2.sli.llf+dummy1.sli.typ)
summary(reg8)
e<-round(lm.beta(reg8),2)

#adding asterisks to the p-vaules of coefficients of Step 2 regression results
p <- summary(reg8)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
e<-paste(e, mystars, sep="")

#Step1 regression for Inhibition verbal
reg9<-lm(sqrt(vimiVtot.err)~age.test+bas.mat.t+bas.v.std)
summary(reg9)

#Step2 regression for Inhibition verbal
reg10<-lm(sqrt(vimiVtot.err)~age.test+bas.mat.t+bas.v.std+dummy2.sli.llf+dummy1.sli.typ)
summary(reg10)
f<-round(lm.beta(reg10),2)

#adding asterisks to the p-vaules of coefficients of Step 2 regression results
p <- summary(reg10)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
f<-paste(f, mystars, sep="")

#Step1 regression for Inhibition nonverbal
reg11<-lm(sqrt(vimiMtot.err)~age.test+bas.mat.t+bas.v.std)
summary(reg11)

#Step2 regression for Inhibition nonverbal
reg12<-lm(sqrt(vimiMtot.err)~age.test+bas.mat.t+bas.v.std+dummy2.sli.llf+dummy1.sli.typ)
summary(reg12)

library(car)
outlierTest(reg11)#63
outlierTest(reg12)#26
Ureg11<-update(reg11,subset=rownames(a)!="63")
Ureg12<-update(reg12,subset=rownames(a)!="26")
g<-round(lm.beta(reg12),2)

#adding asterisks to the p-vaules of coefficients of Step 2 regression results
p <- summary(Ureg12)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
g<-paste(g, mystars, sep="")
g

#Step1 regression for Planning verbal
reg13<-lm(dkef.st.v.raw~age.test+bas.mat.t+bas.v.std)
summary(reg13)

#Step2 regression for Planning verbal
reg14<-lm(dkef.st.v.raw~age.test+bas.mat.t+bas.v.std+dummy2.sli.llf+dummy1.sli.typ)
summary(reg14)

outlierTest(reg13)#77
outlierTest(reg14)#27
plot(reg13,which=c(1,2,5))#77,27,2
plot(reg14,which=c(1,2,5))
Ureg13<-update(reg13,subset=rownames(a)!="77")
Ureg14<-update(reg14,subset=rownames(a)!="27")

h<-round(lm.beta(Ureg14),2)

#adding asterisks to the p-vaules of coefficients of Step 2 regression results
p <- summary(reg14)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
h<-paste(h, mystars, sep="")
h
#Step1 regression for Planning nonverbal
reg15<-lm(dkef.st.nv.raw~age.test+bas.mat.t+bas.v.std)
summary(reg15)

#Step2 regression for Planning nonverbal
reg16<-lm(dkef.st.nv.raw~age.test+bas.mat.t+bas.v.std+dummy2.sli.llf+dummy1.sli.typ)
summary(reg16)
i<-round(lm.beta(reg16),2)

#adding asterisks to the p-vaules of coefficients of Step 2 regression results
p <- summary(reg16)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
i<-paste(i, mystars, sep="")

#Step1 regression for Switching verbal
reg17<-lm(TMswitch.cost~age.test+bas.mat.t+bas.v.std)
summary(reg17)

#Step2 regression for Switching verbal
reg18<-lm(TMswitch.cost~age.test+bas.mat.t+bas.v.std+dummy2.sli.llf+dummy1.sli.typ)
summary(reg18)
j<-round(lm.beta(reg18),2)

#adding asterisks to the p-vaules of coefficients of Step 2 regression results
p <- summary(reg18)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
j<-paste(j, mystars, sep="")

#Step1 regression for Switching nonverbal
reg19<-lm(ied.tot.err~age.test+bas.mat.t+bas.v.std)
summary(reg19)

#Step2 regression for Switching nonverbal
reg20<-lm(ied.tot.err~age.test+bas.mat.t+bas.v.std+dummy2.sli.llf+dummy1.sli.typ)
summary(reg20)
k<-round(lm.beta(reg20),2)

#adding asterisks to the p-vaules of coefficients of Step 2 regression results
p <- summary(reg20)$coefficients[-1,4]  
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
k<-paste(k, mystars, sep="")

#a matrix of the standardised beta values of Step 2 regression for each of 10 variables
Mat<-rbind(b,c,d,e,f,g,h,i,j,k)
Mat

#adding column 1 for "Total R^2 accounted for by the model" & column 7 for "��R^2 Step2"
Mat<-data.frame(a=c(0,0,0,0,0,0,0,0,0,0),Mat,a=c(0,0,0,0,0,0,0,0,0,0))
Mat

#column1: R^2 from summaries of Step 2 regressions
Mat[,1]<-c(round(summary(reg2)$r.squared,2),
round(summary(reg4)$r.squared,2),
round(summary(reg6)$r.squared,2),
round(summary(reg8)$r.squared,2),
round(summary(reg10)$r.squared,2),
round(summary(Ureg12)$r.squared,2),
round(summary(Ureg14)$r.squared,2),
round(summary(reg16)$r.squared,2),
round(summary(reg18)$r.squared,2),
round(summary(reg20)$r.squared,2))
Mat

#column 7: the difference of R squared between Step 1 and Step 2 for each of 10 variables

#anova test to find the difference of R squared is significant
anova(reg1,reg2)

#adding asterisks to the p-value of anova test
p <- anova(reg1,reg2)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p1<-paste(round((summary(reg2)$r.squared)-(summary(reg1)$r.squared),2), mystars, sep="")

#anova test to find the difference of R squared is significant
anova(reg3,reg4)

#adding asterisks to the p-value of anova test
p <- anova(reg3,reg4)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p2<-paste(round((summary(reg4)$r.squared)-(summary(reg3)$r.squared),2), mystars, sep="")

#anova test to find the difference of R squared is significant
anova(reg5,reg6)

#adding asterisks to the p-value of anova test
p <- anova(reg5,reg6)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p3<-paste(round((summary(reg6)$r.squared)-(summary(reg5)$r.squared),2), mystars, sep="")

#anova test to find the difference of R squared is significant
anova(reg7,reg8)
p <- anova(reg7,reg8)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p4<-paste(round((summary(reg8)$r.squared)-(summary(reg7)$r.squared),2), mystars, sep="")

#anova test to find the difference of R squared is significant
anova(reg9,reg10)

#adding asterisks to the p-value of anova test
p <- anova(reg9,reg10)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p5<-paste(round((summary(reg10)$r.squared)-(summary(reg9)$r.squared),2), mystars, sep="")

#anova test to find the difference of R squared is significant
anova(Ureg11,Ureg12)

#adding asterisks to the p-value of anova test
p <- anova(Ureg11,Ureg12)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p6<-paste(round((summary(reg12)$r.squared)-(summary(reg11)$r.squared),2), mystars, sep="")
p6
#anova test to find the difference of R squared is significant
anova(Ureg13,Ureg14)

#adding asterisks to the p-value of anova test
p <- anova(Ureg13,Ureg14)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p7<-paste(round((summary(reg14)$r.squared)-(summary(reg13)$r.squared),2), mystars, sep="")

#anova test to find the difference of R squared is significant
anova(reg15,reg16)

#adding asterisks to the p-value of anova test
p <- anova(reg15,reg16)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p8<-paste(round((summary(reg16)$r.squared)-(summary(reg15)$r.squared),2), mystars, sep="")

#anova test to find the difference of R squared is significant
anova(reg17,reg18)

#adding asterisks to the p-value of anova test
p <- anova(reg17,reg18)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p8<-paste(round((summary(reg16)$r.squared)-(summary(reg15)$r.squared),2), mystars, sep="")

#anova test to find the difference of R squared is significant
anova(reg19,reg20)

#adding asterisks to the p-value of anova test
p <- anova(reg19,reg20)[-1,6] 
mystars <- ifelse(p < 0.001, "***", ifelse(p < 0.01, "** ", ifelse(p < 0.05, "*  ", ifelse(p<0.1,".  ","   ")))) 
p10<-paste(round((summary(reg20)$r.squared)-(summary(reg19)$r.squared),2), mystars, sep="")



#column 7:Step 2 explains more of the overall variance than Step 1 model
#��R^2 = R^2 of Step2 model - R^2 of Step1 model
#a matrix of the difference of R squared between Step 1 and Step 2 for each of 10 variables
Mat[,7]<-c(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10)
Mat

#adding column/row names
colnames(Mat)<-c("Total R^2 accounted for by the model","age","NVIQ","VIQ","SLI-vs.-LLF","SLI-vs.-typical","R^2 step2")
rownames(Mat)<-c("ELWM verbal","ELWM nonverbal","Fluency verbal","Fluency nonverbal","Inhibition verbal","Inhibition nonverbal","Planning verbal","Planning nonverbal","Switching verbal","Switching nonverbal")

Mat





#TABLE4########################################################################################################################################################################################################################
#subset participants from group T in order to extract children who are in the same age range as those in group S
#participants in group S are aged between 8 and 14 years old
#extract children who are aged more than 96 months
ageT<-subset(T,age.test>=96)
stem(ageT$age.test)#shows children who are aged between 8 and 14 years old
dim(ageT)#63 children are chosen


#percentage of children in group S, whose scores were 1 sd below the mean of 10 variables for children in group ageT
b<-c(round(sum(S$wmtbc.lr.raw<mean(ageT$wmtbc.lr.raw)-sd(ageT$wmtbc.lr.raw))/41*100,0),round(sum(S$ooo.total<mean(ageT$ooo.total)-sd(ageT$ooo.total))/41*100,0),
round(sum(S$verbal.fluency<mean(ageT$verbal.fluency)-sd(ageT$verbal.fluency))/41*100,0),round(sum(S$design.fluency<mean(ageT$design.fluency)-sd(ageT$design.fluency))/41*100,0),
round(sum(S$vimiVtot.err>mean(ageT$vimiVtot.err)+sd(ageT$vimiVtot.err))/41*100,0),round(sum(S$vimiMtot.err>=mean(ageT$vimiMtot.err)+sd(ageT$vimiMtot.err))/41*100,0),
round(sum(S$dkef.st.v.raw<mean(ageT$dkef.st.v.raw)-sd(ageT$dkef.st.v.raw))/41*100,0),round(sum(S$dkef.st.nv.raw<mean(ageT$dkef.st.nv.raw)-sd(ageT$dkef.st.nv.raw))/41*100,0)
,round(sum(S$TMswitch.cost>mean(ageT$TMswitch.cost)+sd(ageT$TMswitch.cost))/41*100,0),round(sum(S$ied.tot.err>mean(ageT$ied.tot.err)+sd(ageT$ied.tot.err))/41*100,0))

b

#percentage of children in group S, whose scores were 2 sd below the means of 10 variables for children in group ageT
c<-c(round(sum(S$wmtbc.lr.raw<mean(ageT$wmtbc.lr.raw)-2*sd(ageT$wmtbc.lr.raw))/41*100,0),round(sum(S$ooo.total<mean(ageT$ooo.total)-2*sd(ageT$ooo.total))/41*100,0),
round(sum(S$verbal.fluency<mean(ageT$verbal.fluency)-2*sd(ageT$verbal.fluency))/41*100,0),round(sum(S$design.fluency<mean(ageT$design.fluency)-2*sd(ageT$design.fluency))/41*100,0),
round(sum(S$vimiVtot.err>mean(ageT$vimiVtot.err)+2*sd(ageT$vimiVtot.err))/41*100,0),round(sum(S$vimiMtot.err>mean(ageT$vimiMtot.err)+2*sd(ageT$vimiMtot.err))/41*100,0),
round(sum(S$dkef.st.v.raw<mean(ageT$dkef.st.v.raw)-2*sd(ageT$dkef.st.v.raw))/41*100,0),round(sum(S$dkef.st.nv.raw<mean(ageT$dkef.st.nv.raw)-2*sd(ageT$dkef.st.nv.raw))/41*100,0)
,round(sum(S$TMswitch.cost>mean(ageT$TMswitch.cost)+2*sd(ageT$TMswitch.cost))/41*100,0),round(sum(S$ied.tot.err>mean(ageT$ied.tot.err)+2*sd(ageT$ied.tot.err))/41*100,0))

c

#a matrix of b and c
m<-cbind(b,c)
m

#adding column/row names
colnames(m)<-c("<1SD","<2SD")
rownames(m)<-c("ELWM verbal","ELWM nonverbal","Fluency verbal","Fluency nonverbal","Inhibition verbal","Inhibition nonverbal","Planning verbal","Planning nonverbal","Switching verbal","Switching nonverbal")
m

#Sample size, power, and effect size ###################################################################################################################################################################################

library(pwr)

#required sample size if the model is to have 0.15 effect size and 0.8 power at sig. level of 0.05
pwr.f2.test(u=5, f2=0.15,sig.level=.05,power=.8)

#power when we have 160 participants, 5 predictors and 0.15 effect size at sig.level of 0.05
pwr.f2.test(u=5,v=160, f2=0.15,sig.level=.05)

#power when we have 120 participants, 5 predictors and 0.15 effect size at sig.level of 0.05
pwr.f2.test(u=5,v=120, f2=0.15,sig.level=.05)















